Datos espaciales en R

Omar E. Barrantes Sotela

Escuela de Ciencias Geográficas

2022-11-12

Apertura de tabla de datos

En R es posible cargar una tabla de datos con múltiples variables a estas se les conoce como data.frame(). Es importante conocer algunos detalles de la tabla de datos.

accidentes <- read.csv("./datos/COSEVI_accidentes_complete.csv",
                       fileEncoding="utf8",
                       sep = ",", # Separador de campos
                       dec = ".", # Simbolo decimal
                       header = TRUE, # Tabla tiene encabezados?
                       na.strings = "NA")

Librerias para datos espaciales

En R es posible abrir diversos tipos de datos espaciales. Primero se va a cargar las siguientes librerías:

library(sf)
library(sp)
library(rgdal)
#Otras librerias
library(tidyverse)
library(ggplot2)

Datos puntuales referenciados

Los datos corresponden a una locación específica (p.e., latitud-longitud; x-y). Pueden ser sitios, eventos, etc. En este caso se usará la tabla de datos de accidentes de tránsito en el Cantón de Heredia.

act.sp <- sf::st_as_sf(accidentes,
                       coords = c("NX_WGS84","NY_WGS84"))
st_crs(act.sp) <- 4326 # WGS84 EPSG code 4326

Polígonos

En este caso podemos cargar los datos shapefile con los límites de distrito de Costa Rica, y descargados del SNIT, y emplear solo un subconjunto de estos. En este caso solo interesa los distritos en la GAM de Heredia (por lo que se excluye el distrito Vara Blanca).

dist <- read_sf("./datos/shp/IGN_5_limitedistrital_5k_2021.shp")
heredia <- subset(dist,
                  canton == "Heredia" & codigo_dta < 40105)

Operaciones básicas entre datos espaciales

Es posible realizar operaciones como una intersección entre capas vectoriales.

proj.dat <- sf::st_crs(act.sp)
if (proj.dat$input != "EPSG:8908"){
  act.sp <- st_transform(act.sp,8908)  # Proyecta al CR-SIRGAS
}
act.sp.401 <- st_intersection(act.sp, heredia)
#st_crs(heredia)

Importante

En caso de que las capas presenten diferentes proyecciones, se utiliza la condicional if para verificar esa condición, y proyectar la capa de interés en caso de que no cumpla ese criterio.

Visualización básica

Se crea un mapa básico con los accidentes de tránsito en el Cantón de Heredia y se asigna a una variable para incluso agregar más elementos.

library(ggforce)
library(ggspatial)
library(showtext)

plt01 <- ggplot()+
  annotation_spatial(heredia)+
  layer_spatial(act.sp.401, aes(col= CodAcc2)) +
  labs(colour="Tipo") +
  theme_bw()

Heredia: Mapa de accidentes

plt01

Figura 1: Heredia: Accidentes de tránsito por tipo

Agrupación de datos

Es posible usar una estrategia de agrupación de los accidentes por medio de K Medias, según la variable total de involucrados:

library(cluster)

acc.clus01 <- kmeans(c(act.sp.401$hora),
                     4, # Cantidad de grupos
                     nstart = 20 # Mínimo de elementos en los grupos
                     )
acc.clus01
K-means clustering with 4 clusters of sizes 480, 488, 123, 399

Cluster means:
       [,1]
1 14.758333
2  8.323770
3  2.520325
4 19.889724

Clustering vector:
   [1] 4 2 2 2 1 3 3 3 2 1 2 2 1 2 3 2 4 4 4 2 1 2 4 1 2 3 4 4 3 3 2 2 1 2 3 2 1
  [38] 2 4 3 1 2 3 4 2 4 3 4 2 1 2 4 4 1 2 1 3 2 4 1 2 4 1 2 2 2 2 1 2 1 4 1 1 2
  [75] 1 2 4 2 1 2 1 2 4 1 4 4 4 2 1 1 1 4 1 2 2 2 2 4 1 2 2 1 4 4 1 1 1 4 2 2 2
 [112] 4 2 4 4 1 2 4 2 4 2 4 2 2 1 2 2 2 1 2 2 2 2 3 1 2 4 4 2 2 1 2 1 1 4 1 4 1
 [149] 2 4 1 2 1 4 1 2 1 4 2 4 4 2 4 2 1 2 1 4 2 2 1 4 2 1 4 1 2 4 4 2 2 2 2 2 2
 [186] 2 4 1 2 1 2 4 1 4 2 2 2 3 2 2 1 1 2 1 1 4 2 4 2 2 1 2 2 2 1 4 2 2 4 2 4 1
 [223] 1 2 2 4 1 2 2 3 4 4 2 4 4 1 4 2 1 3 3 1 2 1 4 2 3 2 2 4 4 2 4 2 2 1 4 1 4
 [260] 1 1 4 4 4 1 4 2 4 4 2 2 4 4 1 2 3 1 4 1 1 4 2 1 4 2 2 1 4 2 3 2 1 1 2 4 2
 [297] 4 1 2 4 2 2 2 2 4 1 2 3 1 2 2 1 2 1 2 2 2 1 4 4 4 1 2 1 1 3 3 1 1 1 2 2 4
 [334] 2 2 4 3 2 2 1 3 1 1 4 4 3 4 4 1 4 4 2 1 1 4 2 4 4 4 1 1 4 1 1 2 1 4 1 4 1
 [371] 4 1 1 4 4 2 1 4 1 2 4 2 1 1 2 2 2 2 2 2 2 1 2 4 1 1 1 3 3 4 4 2 2 2 1 3 3
 [408] 1 2 3 2 2 4 1 3 1 2 1 4 2 1 3 3 1 1 4 4 4 2 2 2 2 2 4 1 2 1 3 4 1 4 4 1 1
 [445] 4 1 1 4 1 3 3 1 2 3 4 3 1 4 3 4 2 4 1 2 4 4 2 4 2 2 1 1 1 4 2 2 1 1 4 4 1
 [482] 3 4 2 2 1 4 1 2 1 1 4 1 2 4 1 2 2 2 4 2 2 2 1 1 4 2 4 4 4 3 1 4 2 3 3 4 4
 [519] 1 4 2 2 4 2 2 4 1 4 2 4 2 1 1 1 1 4 4 2 2 1 4 4 1 2 1 2 2 3 2 1 4 4 2 2 1
 [556] 1 2 4 4 3 3 4 2 2 2 4 2 4 2 1 4 2 4 1 4 1 1 2 2 2 1 2 4 4 4 4 3 1 2 1 1 4
 [593] 4 1 4 2 4 1 4 4 2 3 2 2 4 3 1 2 1 4 4 2 1 2 1 4 2 1 1 2 2 2 3 3 4 1 4 2 1
 [630] 4 4 2 2 1 2 1 4 3 1 4 1 1 3 2 2 4 2 1 2 4 4 3 4 4 1 4 4 2 2 2 4 4 4 4 1 3
 [667] 1 3 4 3 4 4 4 1 1 4 2 4 4 4 2 1 2 3 2 1 2 4 2 2 4 1 1 2 1 4 1 4 4 2 4 1 4
 [704] 1 4 2 1 4 4 1 4 4 1 2 4 3 2 1 4 2 3 1 2 3 3 2 4 1 4 1 4 4 1 4 4 4 2 2 1 1
 [741] 4 2 2 2 2 2 1 4 2 4 4 1 1 1 1 2 4 2 3 3 4 1 1 1 2 2 1 2 2 4 2 2 1 4 4 4 4
 [778] 1 4 1 4 4 1 2 2 1 2 3 3 1 4 1 1 1 1 2 1 1 2 1 1 1 1 2 4 2 1 4 3 3 3 4 3 3
 [815] 4 4 4 4 2 4 2 2 2 4 1 1 2 1 2 1 4 1 1 2 3 1 4 1 1 1 2 2 1 2 2 1 1 1 4 4 2
 [852] 2 1 4 2 2 1 1 1 2 2 4 1 2 1 1 4 1 1 4 3 3 3 4 2 1 2 2 4 4 2 1 4 3 4 3 4 1
 [889] 4 2 1 1 2 1 2 4 1 3 3 4 3 2 1 4 2 4 1 1 2 1 1 1 1 3 2 1 1 3 4 4 1 1 4 1 2
 [926] 2 1 2 1 4 1 4 1 1 2 1 1 1 2 3 1 2 1 1 2 1 1 4 2 2 2 2 2 2 1 1 2 1 4 1 2 1
 [963] 1 1 1 1 2 4 1 1 1 2 1 2 3 2 1 1 1 2 1 1 4 1 1 4 2 1 1 4 1 2 4 3 4 3 4 4 4
[1000] 2 1 2 2 2 2 2 2 4 1 1 4 2 2 1 1 3 4 2 4 2 4 1 4 1 1 4 2 1 4 2 1 4 2 1 1 1
[1037] 1 4 1 3 1 4 1 1 1 2 1 2 2 1 4 1 2 1 1 1 1 2 1 1 3 2 2 1 1 4 4 2 2 2 2 1 2
[1074] 1 2 2 1 1 2 2 1 1 2 2 1 1 1 3 2 2 1 2 1 2 1 2 2 4 2 1 1 4 1 1 1 2 2 4 1 1
[1111] 4 4 2 2 1 1 2 2 2 1 3 4 4 4 4 2 3 1 1 1 1 4 2 2 4 1 4 4 4 1 1 1 2 1 1 1 1
[1148] 4 2 1 2 2 4 2 1 2 4 1 1 1 1 1 4 2 1 2 2 1 4 2 2 2 2 2 2 1 1 1 2 4 1 4 4 4
[1185] 4 1 3 4 2 2 2 2 2 1 1 4 1 2 3 3 3 3 2 2 3 3 2 3 3 4 4 2 3 3 1 3 4 2 1 2 1
[1222] 2 1 2 2 1 2 1 1 1 2 4 1 1 2 1 4 2 4 2 2 1 1 1 2 1 2 2 1 4 2 1 1 2 2 1 2 4
[1259] 1 4 1 1 2 4 1 4 1 2 1 4 4 1 4 1 2 1 2 2 4 1 1 1 2 1 1 4 4 1 4 2 1 2 1 4 1
[1296] 2 2 1 1 2 3 2 3 1 4 2 2 4 2 4 2 2 2 4 2 4 2 2 4 2 2 4 2 4 4 1 4 1 2 3 4 3
[1333] 2 1 1 1 4 2 1 1 1 2 4 4 2 1 3 1 1 2 4 3 1 2 4 3 1 4 1 1 4 3 2 2 2 2 1 1 3
[1370] 2 1 2 2 2 2 1 4 4 2 4 4 2 1 1 4 1 3 4 2 4 2 2 1 4 4 4 3 1 2 3 2 2 4 1 2 4
[1407] 2 4 2 4 4 2 1 2 2 1 1 3 1 4 4 1 4 4 1 4 1 4 4 2 2 4 2 2 4 1 3 2 2 4 1 4 1
[1444] 1 2 1 4 1 4 1 4 1 2 1 4 1 1 4 2 3 3 4 4 1 4 3 4 3 2 2 1 2 2 4 2 2 4 1 4 1
[1481] 4 1 4 4 2 4 1 2 2 1

Within cluster sum of squares by cluster:
[1] 1267.9667 1324.8443  480.6992 1209.1479
 (between_SS / total_SS =  91.2 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
[6] "betweenss"    "size"         "iter"         "ifault"

Proceso de clasificación

Es posible realizar una visualización de la agrupación, junto con otras variables. En este caso Hora del accidente con respecto al Total de Involucrados.

clus <- as.data.frame(st_drop_geometry(act.sp.401[,c("hora","Total_Involucrados")]))
clusplot(clus[,c("hora","Total_Involucrados")], 
         acc.clus01$cluster,
         color = TRUE,
         shade = TRUE,
         labels = 2,
         lines = 0,
         main = "")

Heredia: Agrupación de accidentes por Hora y Total de Involucrados

Clasificación en la tabla de datos

El resultado de un proceso de agrupación posibilita la elaboración de cartografía temática. En este caso se construyen niveles en la categoría.

act.sp.401$cat_hora <- as.factor(acc.clus01$cluster)
levels(act.sp.401$cat_hora) <- c("Madrugada",
                                 "Matutino",
                                 "Vespertino",
                                 "Nocturno")

Visualización del resultado de la agrupación.

plt02 <- ggplot()+
  annotation_spatial(heredia)+
  layer_spatial(act.sp.401, aes(col= cat_hora)) +
  labs(colour="Tipo Hora") +
  theme_bw()

Heredia: Mapa de accidentes por horario

plt02

Figura 2: Heredia: Accidentes de tránsito agrupado por horario

Selección por atributo

Es posible solo obtener los datos de accidentes agrupados como nocturno y generar una tabla.

library(knitr)
act.sp.401.noct <-  act.sp.401 %>% 
  filter(cat_hora == "Nocturno")
kable(head(act.sp.401.noct[,c("CodAcc2","hora","cat_hora")]))
Tabla 1: Accidentes Nocturnos
CodAcc2 hora cat_hora geometry
Colisión entre vehículos 21 Nocturno POINT (484447.3 1103424)
Colisión entre vehículos 18 Nocturno POINT (485512.7 1103931)
Colisión entre vehículos 19 Nocturno POINT (483727 1102840)
Colisión entre vehículos 18 Nocturno POINT (482787.3 1103283)
Colisión entre vehículos 21 Nocturno POINT (483592.8 1104124)
Colisión entre vehículos 23 Nocturno POINT (483354.4 1102995)

Escritura de datos

Condición que revisa si el archivo existe, caso contrario se escribe en disco.

fileout00 <-  paste0('./datos/shp/actsp401.shp')

if (!file.exists(fileout00)){
  # save as shapefile
  st_write(act.sp.401, 
         dsn = fileout00, 
         driver = 'ESRI Shapefile',
         delete_layer= TRUE)
} else {
  cat("El archivo ya existe! \n")
}
El archivo ya existe!

Muchas gracias

Escuela de Ciencias Geográficas, UNA.